#ifndef _CU_IMAGE_INL_
#define _CU_IMAGE_INL_

//////////////////////////////////////////////////////////////////////////
/// Filtering kernels
//////////////////////////////////////////////////////////////////////////

#define EXP(x)                                                 __expf((x))
#define SIN(x)                                                 __sinf((x))
#define _CLR_DIST(x, y)                                    _ssd3((x), (y))
#define CLR_DIST(x, y)                                      ssd3((x), (y))


texture<uchar4, 2, cudaReadModeNormalizedFloat> g_BoxBlur_inl;

template <int radius> 
__global__ void BoxBlur_kernel(uchar4 * pDst, float inv_area, uint32 w, uint32 h, uint32 p)
{
	int tidx = threadIdx.x + blockIdx.x * blockDim.x;
	int tidy = threadIdx.y + blockIdx.y * blockDim.y;

	if (tidx < w && tidy < h)
	{
		float4 r = {0.0f, 0.0f, 0.0f, 0.0f}; 
		
		for (int ir = -radius; ir <= radius; ir++)
			for (int ic = -radius; ic <= radius; ic++)
			{
				r += tex2D(g_BoxBlur_inl, tidx + 0.5f + ic, tidy + 0.5f + ir);				
			}

			r *= inv_area;

		pDst[tidx + tidy * p] = uc4(r * 255.0f);
	}
}

texture<uchar4, 2, cudaReadModeNormalizedFloat> g_NonLocalMeans_f4_inl;

template <int radius, int window>
__global__ void NonLocalMeans_tex_kernel(uchar4 * pDst, 
																		 float    inv_sigma_sq, 
																		 float    inv_area, 
																		 float    inv_window_area, 
																		 float    inv_noise, 
																		 float    weight_threshold, 
																		 uint32 w, uint32 h, uint32 p)
{
	const int tidx = threadIdx.x + blockIdx.x * blockDim.x;
	const int tidy = threadIdx.y + blockIdx.y * blockDim.y;

	float tx = tidx + 0.5f;
	float ty = tidy + 0.5f;

	if (tidx < w && tidy < h)
	{
		float4 r = {0.0f, 0.0f, 0.0f, 0.0f};

		float  weight_sum = 0.0f;
		float  weight = 0.0f;

		float  weight_threshold_counter = 0.0f;
		
		for (float ir = -radius; ir <= radius; ir++)
			for (float ic = -radius; ic <= radius; ic++)
			{							
				weight = 0.0f;
				
				#pragma unroll
				for (float j = - window; j <= window; j++)
					#pragma unroll
					for (float i = -window; i <= window; i++)
					{
						float4 c_r_c = tex2D(g_NonLocalMeans_f4_inl, tx + i+0,  ty + j+0);
						float4 c_j_i = tex2D(g_NonLocalMeans_f4_inl, tx + i+ic, ty + j+ir);

						weight += CLR_DIST(c_r_c, c_j_i);
					}				

				weight = weight * inv_window_area;

				weight = EXP( -(weight * inv_noise) - ((ir*ir)+(ic*ic)) * inv_sigma_sq ); // Gaussian component
				
				weight_threshold_counter += (weight >= weight_threshold ? 1.0f : 0.0f);

				weight_sum += weight;

				r += tex2D(g_NonLocalMeans_f4_inl, tx + ic, ty + ir) * weight;
			}

		r /= weight_sum;

		float4 c00 = tex2D(g_NonLocalMeans_f4_inl, tx, ty);		
		r = lerp(c00, r, weight_threshold_counter * inv_area);

		pDst[tidx + tidy * p] = uc4(r * 255.0f);
	}
}

SD_CALL  
uint _sad3(const uchar4 & u, const uchar4 & v) 
{ 
	return __usad(u.x, v.x, __usad(u.y, v.y, __usad(u.z, v.z, 0))); // Sum of Absolute Differences for xyz components only using intrinsic
} 

texture<uchar4, 2, cudaReadModeElementType> g_NonLocalMeans_uc4_inl;

extern __shared__ uchar4 smem_nonlocalmeans_inl[];

template <int radius, int window>
__global__ void NonLocalMeans_tex_smem_kernel(uchar4 * pDst, 
																						  float    inv_sigma_sq, 
																							float    inv_area, 
																							float    inv_window_area, 
																							float    inv_noise, 
																							float    weight_threshold, 
																							uint32 xdim, uint32 ydim,
																							uint32 r, uint32 c, 																						
																							uint32 w, uint32 h, uint32 p)
{
	int tidx = threadIdx.x + blockIdx.x * blockDim.x;
	int tidy = threadIdx.y + blockIdx.y * blockDim.y;	

	int xbase = threadIdx.x;
	int ybase = threadIdx.y;

	int xcoord = tidx - window - radius;
	int ycoord = tidy - window - radius;

	for (int ir = 0; ir < r; ir++)
	{
		xbase = threadIdx.x;
		xcoord = tidx - window - radius;

		for (int ic = 0; ic < c; ic++)
		{
			if (xbase < xdim && ybase < ydim)
			{
				smem_nonlocalmeans_inl[xbase + __mul24(ybase, xdim)] = tex2D(g_NonLocalMeans_uc4_inl, xcoord+0.5f, ycoord+0.5f);
			}

			xbase += blockDim.x;
			xcoord += blockDim.x;	
		}
		ybase += blockDim.y;
		ycoord += blockDim.y;
	}

	xbase = radius + window + threadIdx.x;
	ybase = radius + window + threadIdx.y;

	__syncthreads();

	if (tidx < w && tidy < h)
	{
		float4 r = {0.0f, 0.0f, 0.0f, 0.0f};

		float  weight_sum = 0.0f;
		float  weight = 0.0f;

		const float rcp_255 = 1.0f / 255.0f;
		const float rcp_255_sq = rcp_255 * rcp_255;
		float weight_threshold_counter = 0.0f;
	
		for (int ir = -radius; ir <= radius; ir++)
			for (int ic = -radius; ic <= radius; ic++)
			{							
				weight = 0.0f;

				#pragma unroll
				for (int j = - window; j <= window; j++)
				{					
					int j1 = __mul24((ybase+j+ 0), xdim);
					int j2 = __mul24((ybase+j+ir), xdim);
					#pragma unroll	
					for (int i = -window; i <= window; i++)
					{
						float4 c_r_c = f4(smem_nonlocalmeans_inl[xbase+i+ 0 + j1]);
						float4 c_j_i = f4(smem_nonlocalmeans_inl[xbase+i+ic + j2]); 

						weight += CLR_DIST(c_r_c, c_j_i);
					}				
				}

				weight = weight * inv_window_area * rcp_255_sq;

				weight = EXP( -(weight * inv_noise) - (__mul24(ir,ir)+__mul24(ic,ic)) * inv_sigma_sq ); // Gaussian component
				
				weight_threshold_counter += (weight >= weight_threshold ? 1.0f : 0.0f);

				weight_sum += weight;

				r += smem_nonlocalmeans_inl[xbase+ic+ __mul24((ybase+ir), xdim)] * weight * rcp_255;
			}

		r /= weight_sum;

		float4 c00 = smem_nonlocalmeans_inl[xbase + __mul24((ybase), xdim)] * rcp_255;
		r = lerp(c00, r, weight_threshold_counter * inv_area);

		pDst[tidx + __mul24(tidy, p)] = uc4(r * 255.0f);
	}
}

template <int radius, int window>
__global__ void NonLocalMeans_gmem_smem_kernel(uchar4 * pDst, 
																							 uchar4 * pSrc,
																						   float    inv_sigma_sq, 
																							 float    inv_area, 
																							 float    inv_window_area, 
																							 float    inv_noise, 
																							 float    weight_threshold, 
																							 int32 xdim, int32 ydim,
																							 int32 r, int32 c, 																						
																							 int32 w, int32 h, int32 p)
{
	int tidx = threadIdx.x + blockIdx.x * blockDim.x;
	int tidy = threadIdx.y + blockIdx.y * blockDim.y;	

	int xbase = threadIdx.x;
	int ybase = threadIdx.y;

	int xcoord = tidx - window - radius;
	int ycoord = tidy - window - radius;

	for (int ir = 0; ir < r; ir++)
	{
		xbase = threadIdx.x;
		xcoord = tidx - window - radius;

		for (int ic = 0; ic < c; ic++)
		{
			if ((xbase < xdim) && (ybase < ydim))
			{
				smem_nonlocalmeans_inl[xbase + __mul24(ybase, xdim)] = pSrc[ min(max(xcoord, 0), w-1) + __mul24(min(max(ycoord, 0), h-1), p) ];
			}

			xbase += blockDim.x;
			xcoord += blockDim.x;	
		}
		ybase += blockDim.y;
		ycoord += blockDim.y;
	}

	xbase = radius + window + threadIdx.x;
	ybase = radius + window + threadIdx.y;

	__syncthreads();

	if (tidx < w && tidy < h)
	{
		float4 r = {0.0f, 0.0f, 0.0f, 0.0f};

		float  weight_sum = 0.0f;
		float  weight = 0.0f;

		float rcp_255 = 1.0f / 255.0f;
		float rcp_255_sq = rcp_255 * rcp_255;

		float weight_threshold_counter = 0.0f;
	
		for (int ir = -radius; ir <= radius; ir++)
			for (int ic = -radius; ic <= radius; ic++)
			{							
				weight = 0.0f;

				#pragma unroll
				for (int j = - window; j <= window; j++)
				{
					int j1 = __mul24((ybase+j+ 0), xdim);
					int j2 = __mul24((ybase+j+ir), xdim);
					#pragma unroll	
					for (int i = -window; i <= window; i++)
					{
						float4 c_r_c = f4(smem_nonlocalmeans_inl[xbase+i+ 0 + j1]);
						float4 c_j_i = f4(smem_nonlocalmeans_inl[xbase+i+ic + j2]); 

						weight += CLR_DIST(c_r_c, c_j_i);
					}				
				}

				weight = weight * inv_window_area * rcp_255_sq;

				weight = EXP( -(weight * inv_noise) - (__mul24(ir,ir)+__mul24(ic,ic)) * inv_sigma_sq ); // Gaussian component
				
				weight_threshold_counter += (weight >= weight_threshold ? 1.0f : 0.0f);

				weight_sum += weight;

				r += smem_nonlocalmeans_inl[xbase+ic+ __mul24((ybase+ir), xdim)] * weight * rcp_255;
			}

		r /= weight_sum;

		float4 c00 = smem_nonlocalmeans_inl[xbase + __mul24(ybase, xdim)] * rcp_255;
		r = lerp(c00, r, weight_threshold_counter * inv_area);

		pDst[tidx + tidy * p] = uc4(r * 255.0f);
	}
}

texture<uchar4, 2, cudaReadModeNormalizedFloat> g_Lanczos_inl;

template <int r>
__device__ float lanczos(float x)
{
	const float m_pi = 3.14159265358979f;

	float result = 0.0f;

  if (x >= -r && x <= r) 
	{  
		float a = x * m_pi; // PI value
		float b = (r * (SIN(a / r)) * (SIN(a)) / (a * a));

		result = (x == 0.0f) ? 1.0f : b;
	}

  return result;
}

template <int radius>
__global__ void Lanczos_kernel(uchar4 *pDst, float factor, float blur, float support, float scale, uint32 w, uint32 h, uint32 p)
{
	int tidx = threadIdx.x + blockDim.x * blockIdx.x;
	int tidy = threadIdx.y + blockDim.y * blockIdx.y;

	if(tidx < w && tidy < h)
	{
		float weight_sum = 0.0f;
		float weight = 0.0f;
		float4 r = { 0.0f, 0.0f, 0.0f, 0.0f };

    float center = tidx / factor;
    int32 start = (int32) max(center-support+0.5f, (float)0);
    int32 stop  = (int32) min(center+support+0.5f, (float)w);
    
		int32 nmax = stop-start;
    
		float s = start - center;

    for (int n=0; n<nmax; ++n, ++s)
    {
      weight = lanczos<radius> (s * scale);
      weight_sum += weight;
			r += (tex2D(g_Lanczos_inl, tidy + 0.5f, start+n + 0.5f)*weight);
    }

    if (weight_sum != 0.0f)
    {
      /* Normalize. */
      r /= weight_sum;
    }

#if RESAMPLING_ADD_CONTRAST
		r = pow(r, RESAMPLING_POWER);
#endif

		pDst[tidx + tidy * p] = uc4(clamp(r*255.0f, 0, 255.0f));
	}
}

//////////////////////////////////////////////////////////////////////////////
/////// Filtering wrappers
//////////////////////////////////////////////////////////////////////////////

template <int radius>
void BoxBlur(cuImage & dst, cuImage & src)
{
	cudaError_t err;

	uint32 w = src.Width();
	uint32 h = src.Height();
	uint32 p = src.Pitch();

	err = cudaBindTextureToArray(g_BoxBlur_inl, src.m_data.aPtr());

	if (err != cudaSuccess) std::cerr << "Error binding texture to array in BoxBlur<R>" << std::endl;

	dim3 block(32, 8);
	dim3 grid( w / block.x + ((w % block.x) ? 1: 0), 
		         h / block.y + ((h % block.y) ? 1: 0) );

	float inv_area = 1.0f / sqr(2.0f*radius+1.0f);

	BoxBlur_kernel<radius><<<grid, block>>>(dst.m_data.dPtr(), 
		                                      inv_area, 
																					w, h, p);

	err = cudaThreadSynchronize();

	if (err != cudaSuccess) std::cerr << "Error during BoxBlur_kernel<R> execution" << std::endl;
}

template <int radius, int window>
void NonLocalMeans_tex(cuImage & dst, cuImage & src, float noise, float weight_threshold)
{
	cudaError_t err;

	uint32 w = src.Width();
	uint32 h = src.Height();
	uint32 p = src.Pitch();

	err = cudaBindTextureToArray(g_NonLocalMeans_f4_inl, src.m_data.aPtr());

	if (err != cudaSuccess) std::cerr << "Error binding texture to array in NonLocalMeans<R,W>" << std::endl;

	dim3 block(32, 8);
	dim3 grid( w / block.x + ((w % block.x) ? 1: 0), 
		         h / block.y + ((h % block.y) ? 1: 0) );

	float inv_area        = 1.0f / sqr(2.0f*radius+1.0f);
	float inv_sigma_sq    = 1.0f / (2.0f*radius*radius);
	float inv_noise       = 1.0f / noise;
	float inv_window_area = 1.0f / sqr(2.0f*window+1.0f);

	NonLocalMeans_tex_kernel<radius, window><<<grid, block>>>(dst.m_data.dPtr(), 
		                                                    inv_sigma_sq, 
																												inv_area, 
																												inv_window_area, 
																												inv_noise, 
																												weight_threshold, 
																												w, h, p);

	err = cudaThreadSynchronize();

	if (err != cudaSuccess) std::cerr << "Error during NonLocalMeans_kernel<R,W> execution" << std::endl;
}

template <int radius, int window>
void NonLocalMeans_tex_smem(cuImage & dst, cuImage & src, float noise, float weight_threshold)
{
	cudaError_t err;

	uint32 w = src.Width();
	uint32 h = src.Height();
	uint32 p = src.Pitch();

	err = cudaBindTextureToArray(g_NonLocalMeans_uc4_inl, src.m_data.aPtr());
	if (err != cudaSuccess) std::cerr << "Error binding texture to array in NonLocalMeans_tex_smem<R,W>" << std::endl;

	err = cudaBindTextureToArray(g_NonLocalMeans_f4_inl, src.m_data.aPtr());
	if (err != cudaSuccess) std::cerr << "Error binding texture to array in NonLocalMeans_tex_smem<R,W>" << std::endl;

	dim3 block(16, 16);
	dim3 grid( w / block.x + ((w % block.x) ? 1: 0), 
		         h / block.y + ((h % block.y) ? 1: 0) );

	float inv_area        = 1.0f / sqr(2*radius+1);
	float inv_sigma_sq    = 1.0f / (2*radius*radius);
	float inv_noise       = 1.0f / (noise);
	float inv_window_area = 1.0f / sqr(2*window+1);

	int xdim = block.x + 2*window+2*radius;
	int ydim = block.y + 2*window+2*radius;

	int smem = xdim * ydim * sizeof(uchar4);

	int c = xdim / block.x + ((xdim % block.x) ? 1: 0);
	int r = ydim / block.y + ((ydim % block.y) ? 1: 0);

	NonLocalMeans_tex_smem_kernel<radius, window><<<grid, block, smem>>>(dst.m_data.dPtr(), 
		                                                                   inv_sigma_sq, 
																																			 inv_area, 
																																			 inv_window_area, 
																																			 inv_noise, 
																																			 weight_threshold, 
																																			 xdim, ydim, 
																																			 r, c, 
																																			 w, h, p);

	err = cudaThreadSynchronize();

	if (err != cudaSuccess) std::cerr << "Error during NonLocalMeans_tex_smem<R,W> execution" << std::endl;
}

template <int radius, int window>
void NonLocalMeans_gmem_smem(cuImage & dst, cuImage & src, float noise, float weight_threshold)
{
	cudaError_t err;

	uint32 w = src.Width();
	uint32 h = src.Height();
	uint32 p = src.Pitch();

	dim3 block(16, 16);
	dim3 grid( w / block.x + ((w % block.x) ? 1: 0), 
		         h / block.y + ((h % block.y) ? 1: 0) );

	float inv_area        = 1.0f / sqr(2*radius+1);
	float inv_sigma_sq    = 1.0f / (2*radius*radius);
	float inv_noise       = 1.0f / (noise);
	float inv_window_area = 1.0f / sqr(2*window+1);

	int xdim = block.x + 2*window+2*radius;
	int ydim = block.y + 2*window+2*radius;

	int smem = xdim * ydim * sizeof(uchar4);

	int c = xdim / block.x + ((xdim % block.x) ? 1: 0);
	int r = ydim / block.y + ((ydim % block.y) ? 1: 0);

	NonLocalMeans_gmem_smem_kernel<radius, window><<<grid, block, smem>>>(dst.m_data.dPtr(),
		                                                                    src.m_data.dPtr(),
		                                                                    inv_sigma_sq, 
																																			  inv_area, 
																																			  inv_window_area, 
																																			  inv_noise, 
																																			  weight_threshold, 
																																			  xdim, ydim, 
																																			  r, c, 
																																			  w, h, p);

	err = cudaThreadSynchronize();

	if (err != cudaSuccess) std::cerr << "Error during NonLocalMeans_gmem_smem<R,W> execution" << std::endl;
}

template <int radius>
void Lanczos(cuImage & dst, cuImage & src, cuImage & dst_y, float blur)
{
	cudaError_t err;

	uint32 src_w   = src.Width();
	uint32 src_h   = src.Height();
	uint32 src_p   = src.Pitch();

	uint32 dst_w   = dst.Width();
	uint32 dst_h   = dst.Height();
	uint32 dst_p   = dst.Pitch();

	uint32 dst_y_w = dst_y.Width();
	uint32 dst_y_h = dst_y.Height();
	uint32 dst_y_p = dst_y.Pitch();

	float factor_w = dst.Width() / (float) src.Width();
	float factor_h = dst.Height() / (float) src.Height();

	float scale_w = min(factor_w, 1.0f) / blur;
	float scale_h = min(factor_h, 1.0f) / blur;

	float support_w = radius / scale_w;
	float support_h = radius / scale_h;

	if (support_w <= 0.5f) { support_w = 0.5f + 1E-12f; scale_w = 1.0f; }
	if (support_h <= 0.5f) { support_h = 0.5f + 1E-12f; scale_h = 1.0f; }

	err = cudaBindTextureToArray(g_Lanczos_inl, src.m_data.aPtr());

	if (err != cudaSuccess) std::cerr << "Error binding texture to array in Lanczos<R>" << std::endl;

	dim3 block_h(32, 8);
	dim3 grid_h( dst_y_w / block_h.x + ((dst_y_w % block_h.x) ? 1: 0), 
		           dst_y_h / block_h.y + ((dst_y_h % block_h.y) ? 1: 0) );

	Lanczos_kernel<radius><<<grid_h, block_h>>>((uchar4 *) dst_y.m_data.dPtr(), factor_h, blur, support_h, scale_h, dst_y_w, dst_y_h, dst_y_p);

	err = cudaThreadSynchronize();

	if (err != cudaSuccess) std::cerr << "Error during Lanczos_kernel<R> execution (scale h)" << std::endl;

	dst_y.Memcpy(BUFFER_D_TO_A);

	err = cudaBindTextureToArray(g_Lanczos_inl, dst_y.m_data.aPtr());

	if (err != cudaSuccess) std::cerr << "Error binding texture to array in Lanczos<R>" << std::endl;

	dim3 block(32, 8);
	dim3 grid( dst_w / block.x + ((dst_w % block.x) ? 1: 0), 
		         dst_h / block.y + ((dst_h % block.y) ? 1: 0) );

	Lanczos_kernel<radius><<<grid, block>>>((uchar4 *) dst.m_data.dPtr(), factor_w, blur, support_w, scale_w, dst_w, dst_h, dst_p);

	err = cudaThreadSynchronize();

	if (err != cudaSuccess) std::cerr << "Error during Lanczos_kernel<R> execution (scale w)" << std::endl;
}

////////////////////////////////////////////////////////////////////////////
///// Explicit Instantiation
////////////////////////////////////////////////////////////////////////////

template void BoxBlur<5>(cuImage & dst, cuImage & src);
template void NonLocalMeans_tex<3, 3>(cuImage & dst, cuImage & src, float noise, float weight_threshold);
template void NonLocalMeans_tex_smem<3, 3>(cuImage & dst, cuImage & src, float noise, float weight_threshold);
template void NonLocalMeans_gmem_smem<3, 3>(cuImage & dst, cuImage & src, float noise, float weight_threshold);
template void Lanczos<3>(cuImage & dst, cuImage & src, cuImage & dst_y, float blur);

#endif